Probing excitations and cooperatively rearranging regions in deeply supercooled liquids

Upon approaching the glass transition, the relaxation of supercooled liquids is controlled by activated processes, which become dominant at temperatures below the so-called dynamical crossover predicted by Mode Coupling theory (MCT). Two of the main frameworks rationalising this behaviour are dynamic facilitation theory (DF) and the thermodynamic scenario which give equally good descriptions of the available data. Only particle-resolved data from liquids supercooled below the MCT crossover can reveal the microscopic mechanism of relaxation. By employing state-of-the-art GPU simulations and nano-particle resolved colloidal experiments, we identify the elementary units of relaxation in deeply supercooled liquids. Focusing on the excitations of DF and cooperatively rearranging regions (CRRs) implied by the thermodynamic scenario, we find that several predictions of both hold well below the MCT crossover: for the elementary excitations, their density follows a Boltzmann law, and their timescales converge at low temperatures. For CRRs, the decrease in bulk configurational entropy is accompanied by the increase of their fractal dimension. While the timescale of excitations remains microscopic, that of CRRs tracks a timescale associated with dynamic heterogeneity, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${t}^{*} \sim {\tau }_{\alpha }^{0.8}$$\end{document}t*~τα0.8. This timescale separation of excitations and CRRs opens the possibility of accumulation of excitations giving rise to cooperative behaviour leading to CRRs.

Here we use the Roskilde University Molecular Dynamics (RUMD) code. This simulates standard Molecular Dynamics, taking advantage of the many core nature of graphics cards to achieve high performance. In this way, the dynamics simulated are realistic, which is useful for estimating the dynamical quantities of interest here [1,2]. While the dynamics are realistic, the configurations generated are limited to around three decades of increase in relaxation time with respect to the mode-coupling crossover. Unlike, e.g., SWAP which is an unphysical algorithmic improvement to simulate low temperatures, RUMD achieves high performance by optimising the force calculation to utilise thousands of cores available on the GPU. RUMD thus simulates the "true" dynamics of the system as opposed to, e.g., the SWAP algorithm.
We study the binary Lennard-Jones Kob-Andersen (KA) mixture at the non-standard 2:1 and 3:1 composition in the NVT ensemble (Nose-Hoover thermostat) at ρ = 1.400 using the Roskilde University Molecular Dynamics (RUMD) package [3]. We also simulate KA 4:1 for comparing our algorithm's with data in the literature. The interatomic interactions of the KA mixture are defined by with parameters σ AB = 0.80, σ BB = 0.88 and AB = 1.50, BB = 0.50 (α, β = A, B). The pair potential is cut and shifted at r c = 2.5σ αβ . We employ a unit system in which σ AA = 1, AA = 1, and m A = m B = 1 (we refer to this as MD or LJ units). We study system sizes N = 10,002 and 100,002 for KA 2:1 and N = 10,000 and 100,000 for More specifically, the systems were quenched from a high temperature configuration using a thermostat and equilibrated at low temperatures according to the previously mentioned criteria. As the low temperature simulations require absolute stability against round errors occurring over many months it is not recommended to use NVE simulations for the equilibration nor production runs but rather NVT simulations which eliminate any entropic drift in energy due to round off errors.

B. Experimental
Fluorescent dyed-poly methylmethacrylate (PMMA) colloids with a polyhydroxystearic acid comb stabiliser were synthesised using established methods [4]. To enhance spatial resolution between particles, a non-fluorescent PMMA shell was grown on the fluorescent cores and cross-linked with ethylene glycol dimethacrylate (EGDM) to yield a total radius of 270 nm (polydispersity ≈ 8%, Brownian time τ B = 34 ms). PMMA particles were imaged in a density-and refractive index matching mixture of cis-decalin and cyclohexyl bromide.
Samples were loaded into cells constructed of three coverslips on a microscope slide, where two of the coverslips acted as a spacer, and sealed using epoxy. The structural relaxation time was monitored for waiting times of up to 100 days until it reached a steady state at which point the sample was considered to be equilibrated [5]. For example, the state point φ = 0.598 reached equilibrium after 16 days, corresponding to more than 10τ α . Samples were imaged using a Leica SP8 inverted stimulated emission depletion microscope with a 100x oil immersion lens in STED-3D mode, mounted on an optical table. Measurements were taken at least 20 particle diameters from the cell wall to minimise any structural or dynamic influence of surface layering, which typically persisted for around 5 particle diameters. Images were deconvolved with Huygens Professional version 15.05 prior to analysis (Scientific Volume Imaging, The Netherlands, http://svi.nl).

C. Structural Relaxation Time and Dynamic susceptibility
The relaxation time τ α is calculated from the self-part of intermediate scattering function F S (k, t) for a wave vector near the first peak of the A-particle static structure factor. F S (k, t) is fitted with a stretched exponential: where A is a constant and β ≤ 1 is the stretching exponent. Intermediate scattering functions for KA 2:1 and 3:1 are shown in Fig. 1.
The dynamic susceptibility is calculated from the self-part of the overlap function where Here Θ is the Heaviside step function and r a = 0.30. The same value is used for both species of particles. The dynamic susceptibility for KA 2:1 is shown in Fig. 1(c) in the main text.

D. Determining the Mode-Coupling Temperature
The mode-coupling temperature is around T mct = 0.55 for KA 2:1 and 0.7 for KA 3:1 (see Fig. 2). The mode coupling temperature T mct is estimated with a power law fit [6].
where B is a constant. Since the data is expected to go beyond Tmct some data was excluded from the fit as indicated.

E. Excitation Detection Algorithm
To identify the excitations the timescale t a and the length scale a have to be chosen. a has to be small enough that particles can jump the distance but large enough to exclude small fluctuations. After setting a = 0.5 the algorithm is used for different t a . Fig. 3 shows that there is a range where the concentration of excitations does not depend on t a . The smallest t a in this range was chosen (t a = 200 LJ units) as indicated by the black dashed line.
The algorithm is implemented as follows. Trajectories of length t traj = 5t a containing 1000 frames are generated and each particle is tested as a potential excitation. For this the average particle positions of the first and the last t a of the trajectory are compared and if this is smaller than a the particle is rejected. This tests if the particle has committed to a new position (step 1).
If the trajectory has not been rejected we check for each frame if the average positions in the previous t a and the following t a are at least a apart. This makes sure that we exclude particles that move slowly and steadily since we do not consider these to be excitations. This is important at higher temperatures but almost irrelevant at lower temperatures (step 2). To get a first estimate for the hyperbolic tangent fit, for every identified frame f that fulfills the above condition, we go forwards and backwards in time and identify the closest frame for which the distance of the position for the current frame is larger than a/2 (step 3). The time between those two frames is the first guess for ∆t and the central frame of the excitation that produces the shortest δt guess serves as an input for the centre of the fit (step 4). The trajectory is smoothed with a spline (step 5 )and a fitted with a hyperbolic tangent function to extract the exact centre of the excitation and δt (step 6). If this fit fails, the trajectory is rejected as well.
This algorithm becomes less effective when the temperature is raised to around T mct since all particles move considerably on short timescales and excitations are not clearly defined. For low temperatures well below T mct this algorithm is computationally not too expensive since most particles are already rejected in the first step so that the more costly procedure (moving window, fit) are only performed on very few particles.
Excitations where the fit gives a ∆t > t a are excluded since the fit does not make much sense in that case. This algorithm only reliable detects excitations at temperatures lower than T mct since above that excitations (according to our definition) are not as cleanly defined anymore since other types of movement occur (see Fig. 4). Thus we focus on the temperature range T T mct . Since the data in Keys et al. [7] is either out of equilibrium or in the regime T T mct , it is difficult to compare some quantities directly.
Steps in the algorithm for every frame f fullfilling these criteria: δ temp = min(r(t)) r(t)−r(f )>a/2;t − max(r(t)) r(f )−r(t)>a/2;t<f 4. use frame with min δ temp as initial guess for tanh fit 5. smooth trajectory with spline 6. fit trajectory withã tanh ((t − t )/∆t) to get ∆t or reject if fit fails

F. Configurational entropy
The configurational entropy S conf is calculated from the total entropy S by subtracting the basin entropy in the harmonic approximation [8,9] S conf = S − S vib In this equation β = 1/k B T and ω 2 a is an eigenvalue of the Hessian matrix for the inherent state. The average is taken over 200 inherent states with about τ α in between the configurations. Inherent states are found by a simple gradient descent method which works well for the KA liquid.
The total entropy S is calculated via thermodynamic integration from low density (the ideal gas limit) (see Fig. 6). Details of the calculation can be found in Ref. [10]. The definition of the total entropy in MD units requires a value for Planck's constant in MD units. We used h * = 0.1857 which is calculated from using σ = 0.3405·10 −9 m and = 0.03994 kg/mol when converting to MD units [9].
The diffusion coefficient is determined from the long-time limit of the A-particle mean-square displacement. . Yellow shaded regions denotes 5% largest displacements which we use to define particles in CRRs. Jump sizes of particles during an excitation event are indicated as blue crosses. Note that the some of the blue crosses lie at values less that a, the nominal jump displacement. For this there are (at least) two reasons: firstly, a is a fit parameter, the best fit does not always capture the full size of the jump. Secondly it is calculated using the displacement from frame zero, so the part of the trajectory that is fitted could move parallel to the origin so that the fit is much less than a although the displacement between the two is larger than a. 0.08 ± 0.02 0.02 ± 0.01 0.9 ± 0.2 D 21 ± 1 20 ± 1 2.8 ± 0.6 × 10 −3 parabolic fit